c
c    To be linked with wp_xi_cf_hod_12.f
c    This module gives the N1(M),N2(M) and 1-halo term power spectrum
c    designed to model the AGN N1(M) 
c    N1(M) is usual Truncated power-law and assumed to be all satellites.
c    
c    N2(M) is taken from an ascii table gives separate entries for central 
c       and satellite galaxies.
c      
      subroutine init_nm(argc2)      
      ! Read model parameters/ In this model
      ! Mcrit_1,alpha1,filename
      parameter (nmass=50,idv=70)
      character*(*) argc2
      character nm2_fil*80, rline*132, cut_space*132
      real lgrm2(nmass),gn2_cen(nmass),gn2_sat(nmass)
      common/wgt1par/rm_crit1_s,alpha1_s,rm_crit1_c,alpha1_c
      common/wgt2par/nm2,lgrm2,gn2_cen,gn2_sat
c
c     Read Table for population 2
c
      read (argc2,*) rm_crit1_s,alpha1_s,rm_crit1_c,alpha1_c,nm2_fil
      ! For imput, if rm_crit=0 means no such population. In the linear
      ! scale, it is converted to 1.0 so that its log will give 0.
      if (rm_crit1_s.gt.0.and.rm_crit1_s.lt.1e5) then
         rm_crit1_s=10**rm_crit1_s
      else if (rm_crit1_s.le.0) then
         rm_crit1_s=1.0
      end if

      if (rm_crit1_c.gt.0.and.rm_crit1_c.lt.1e5) then
         rm_crit1_c=10**rm_crit1_c
      else if (rm_crit1_c.le.0) then
         rm_crit1_c=1.0
      end if
      !
      open (unit=70,file=nm2_fil,status='old')
      nm2 = 0
      do while (.true.)
         read (70,'(A)',end=3) rline
         rline = cut_space(rline)
         if (index(rline,'!').ne.1) then
            nm2 = nm2+1
            read (rline,*) lgrm2(nm2),gn2_sat(nm2),gn2_cen(nm2)
            if (lgrm2(nm2).gt.1e5) lgrm2(nm2)=log10(lgrm2(nm2))   
!             If entry is lin, covert to log
         end if
      end do        
 3    close (70)
      ! Output Header
      write(6,'(A)')
     $     '! CCF model pl_tbl/Sample1: Truncated PL Sample2: Table'
      write(6,'(A,2(1X,F6.2))')
     $     '! Sample1 (satellite) Log(M_cr1_s) alpha1_s:',
     $     log10(rm_crit1_s),alpha1_s
      write(6,'(A,2(1X,F6.2))')
     $     '! Sample1 (center) Log(M_cr1_c) alpha1_c:',
     $     log10(rm_crit1_c),alpha1_c
      write(6,'(A,A,1X,I2,A)')
     $     '! Sample2: N(M) file:',nm2_fil(1:index(nm2_fil,' ')-1),
     $     nm2,' rows.'
      end
c
c
      real function galn_halo (rm,isam,icen) 
      ! rm: Mass of DMH in solar mass
      ! isam: sample number
      ! icen=0: satellite galaxy, icen=1: central galaxy 
      parameter (nmass=50)
      real lgrm
      real lgrm2(nmass),gn2_cen(nmass),gn2_sat(nmass)
      common/wgt1par/rm_crit1_s,alpha1_s,rm_crit1_c,alpha1_c
      common/wgt2par/nm2,lgrm2,gn2_cen,gn2_sat
      ! By definition, if rm_crit*<=10, no such population.
      if (isam.eq.1.and.icen.eq.0) then ! Sample 1
         if (rm.lt.rm_crit1_s.or.rm_crit1_s.le.10) then ! Satellite
            galn_halo = 0
         else
            galn_halo = (rm/rm_crit1_s)**alpha1_s
         end if
      else if (isam.eq.1.and.icen.eq.1) then
         if (rm.lt.rm_crit1_c.or.rm_crit1_c.le.10) then ! Center
            galn_halo = 0
         else
            galn_halo = (rm/rm_crit1_c)**alpha1_c
         end if
      else ! Sample 2
         lgrm = log10(rm)
         if (lgrm.lt.lgrm2(1)) then
            galn_halo = 0
         else
            imin=1
            imax=nm2
            if (lgrm.gt.lgrm2(nm2)) then ! Above lookup table lim/keep max value
               ff=1.0
               imin = imax-1
            else
               call intra_mono(lgrm2,lgrm,imin,imax,ff)
            end if
            if (icen.eq.0) then ! Satellite
               galn_halo =  (1.-ff)*gn2_sat(imin)+
     $              ff*gn2_sat(imax)
               agaln_halo = max(0.0,galn_halo)
            else ! Center
               galn_halo = (1.-ff)*gn2_cen(imin)+
     $              ff*gn2_cen(imax)
               galn_halo = min(max(0.0,galn_halo),1.0)
            end if
         end if
      end if
      return
      end
c
c
      real function wgt_cf(rm,icen) ! Calculate Weight for
                              ! Expected number of pairs 
                              ! between sample 1 and sample 2
      rm10 = rm/1e10          ! rm in units of 1e10 solar mass/to match with wgt_samp   
      if (icen.eq.0) then     ! Sat vs Sat
         wgt_cf = galn_halo(rm,1,0)*galn_halo(rm,2,0)/(rm10*rm10)
      else   ! Central vs satellite
         wgt_cf = (galn_halo(rm,1,1)*galn_halo(rm,2,0)+
     $        galn_halo(rm,1,0)*galn_halo(rm,2,1))/(rm10*rm10)
      end if   
      return
      end
c
      subroutine model_help()
      print '(A)',' Model Specification: ccf_pl_tbl.f'
      print '(A)',' N(M)_pars: Mcr1_s,alpha1_s,Mcr1_c,alpha1_c,file2'
      print '(A)','  Sample1: N1_s(Mh) =0 (Mh<Mcr1_s); '//
     $     '\\propto(Mh/Mcr1_s)^alpha1_s (>Mcr1_s), satellites'
      print '(A)','  Sample1: N1_c(Mh) =0 (Mh<Mcr1_c); '//
     $     '\\propto(Mh/Mcr1_c)^alpha1_c (>Mcr1_c), center'
      print '(A)','    Mcr=0 means zero N(M) for the population.'
      print '(A)','    Can give Log10(Mcr) instead of Mcr (automatic).'
      print '(A)','  Sample2: ASCII file ''file2''  with 3 cols:'
      print '(A)', '  LogMH [h^{-1}M_sun],N(M)_satellite,N(M)_center'
      print '(A)', ''
      return
      end
